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Using the tight-binding model with long-range Coulomb interactions between electrons, we study 
some of the electronic properties of graphene. The Coulomb interactions are treated with the 
renormalized-ring-diagram approximation. By self-consistently solving the integral equations for 
the Green function, we calculate the spectral density. The obtained result is in agreement with 
experimental observation. In addition, we also compute the density of states, the distribution func- 
tions, and the ground-state energy. Within the present approximation, we find that the imaginary 
part of the self-energy fixed at the Fermi momentum varies as quadratic in energy close to the 
chemical potential, regardless the system is doped or not. This result appears to indicate that the 
electrons in graphene always behave like a moderately correlated Fermi liquid. 

PACS numbers: 71.10.-w, 73.20.At, 81.05.Uw,71.18.+y 
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I. INTRODUCTION 



The graphene is a single-layer honeycomb lattice of 
carbon atoms coated on the surface of some materials^ 2 - 
The Dirac cone structure in the energy spectrum is 
responsible for some of the unusual properties of the 
system^^ Study of the behaviors of electrons in 
graphene is one of the currently focused areas in the 
condensed-matter physics. In the theoretical investiga- 
tions, most of the calculations are based on the continu- 
ous model with simplified Dirac cone dispersion^i and 
the Coulomb interactions between electrons are treated 
in the random-phase approximation (RPA). Since the in- 
teractions are correctly taken into account in the long 
wavelength limit, this approach should reasonably de- 
scribe the low energy behaviors of the electrons within 
the validity of RPA. In such a model, however, the lattice- 
structure effect and the short-range part of the Coulomb 
interaction have been completely neglected. Thus, it is 
desirable to explore this problem by a more realistic ap- 
proach including the effect of the graphene lattice and a 
self-consistent scheme beyond the RPA. 

In this paper, we use the tight- binding model defined 
on the two-dimensional honeycomb lattice to formulate 
the Green function theory of electrons in graphene. In the 
self-energy of electron Green function, the Coulomb in- 
teractions are taken into account with the renormalized- 
ring-diagram approximation (RRDA). This approxima- 
tion is well known to satisfy the microscopic conservation 
laws£ Our recent investigation of the two-dimensional 
electron system- shows that the RRDA can accurately 
reproduce the result of the fixed-node-diffusion Monte 
Carlo simulation for the ground-state energy^ It is 
therefore expected that the RRDA could give more 
reliable description for the behaviors of electrons in 
graphene. 




FIG. 1: (Color online) The structure of a honeycomb lattice. 
The basic vectors of the lattice are ai and a.?. There are 
two sites in each unit cell enclosed by the red lines (solid and 
dotted): black and green. 



II. LATTICE STRUCTURE AND FOURIER 
TRANSFORM 

For readers' convenience, we here briefly review the 
structure of a honeycomb lattice and its reciprocal 
latticei 11 ! 12 For the sake of numerical computation, we 
will also present the mapping between coordinates de- 
fined on the basic vectors of the honeycomb lattice and 
the orthogonal coordinates. 

The graphene lattice is of the honeycomb structure as 
shown in Fig. [TJ A set of basic displacement vectors of 
the lattice is 



ai = (l,0)o 



a 2 



1 y/Z 
4' 2 ; 



(1) 
(2) 



where a is the lattice constant. We will chose a as the 
unit of length, and thereby set a = 1 hereafter. The area 



2 



of the unit cell is 



S 



2 



(3) 



in unit of a 2 = 1. The whole lattice can be viewed as a 
quadrilateral lattice consisting of the unit diamond cells. 
There are two sites in each unit cell: black and green. 
With ai and a 2 , we then define the unit vectors of the 
reciprocal lattice shown in Fig. [2j They are given by 



b! = a 2 x z/S= (l,--j=) 
b 2 = axai/S=(0,4 



(4) 
(5) 



where z is the unit vector in the direction of ai x a 2 . 
The basic displacement vectors of the reciprocal lattice 
are 27rbi and 27rb 2 . The first Brillouin zone (BZ) is the 
hexagon. The diamond enclosed by the red dashed lines 
in Fig. [5] is an equivalent first BZ that is a convenient 
choice for numerical calculation. 

For the use of numerical calculation, we here write 
down the transform between the coordinates on the basis 
of {ai,a2} and the orthogonal axes in real space. Con- 
sider a vector r = (x, y) in the representation of {ai, a 2 }. 
We denote this vector in the orthogonal coordinate sys- 
tem as R = (X, Y) . Then the correspondence between 
these two sets of coordinates is given by 



X = x + | 



Y 



-y- 



(0) 
(7) 



The matrix T of the transform R — Tr is therefore given 
by 



T = 




(8) 



where the first and second columns are the coordinates 
of vectors ai and a 2 respectively. Analogously, in the 
momentum space, we obtain the transform between the 
coordinates of a vector Q defined in the orthogonal sys- 
tem and its projections q on the basis {bi , b 2 }. Denoting 
the transform as Q = Mq, we have 



1 



if 



M = 



1 2 

V3 V3/ 



(9) 



where the first and second columns are respectively the 
vectors bi and b 2 . The two matrices T and M are related 
by M = f'- 1 with T the transpose of f . 

For the later use, we here discuss the Fourier transform. 
The function F(R) defined on the honeycomb lattice sites 
can be expanded as 



F(R) 



BZ 



(10) 
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FIG. 2: (Color online) Reciprocal lattice of honeycomb struc- 
ture. The basic vectors of the lattice are 2nhi and 27rb2. 
The first Brillouin zone is the hexagon with the green bound- 
ary. The diamond enclosed by the red dashed lines is the 
equivalent first Brillouin zone. 



where the Q-integral is over the first BZ with Sbz = 
2(27r) 2 /v / 3 its area, the components of R and Q are 
given in the orthogonal coordinate system, and F(Q) is 
the Fourier component of the function F(R). (A func- 
tion and its Fourier component are distinguished by their 
arguments in this paper. We also adopt the conven- 
tion that a capital vector implies its components defined 
in the orthogonal coordinate system, while a lowercase 
vector means that its components are given in the ba- 
sis {a 1; a 2 } in real space or {b 1; b 2 } in the momentum 
space). In the basis {ai,a 2 }, the function F(R) is given 
by F(Tr) = f{r)- From Eq. (flT)|) , we get the expansion 
for the function /(r), 



BZ 



dq 



F(Mq)e 



tq-r 



(ii) 



Therefore, the Fourier component of / is given by f(q) = 
F(Mq). This relationship is useful for the Fourier trans- 
form of the Coulomb interaction. 



III. TIGHT-BINDING MODEL 

For describing the electron system, we use the tight- 
binding model in which the nearest-neighbor hopping and 
the Coulomb interaction are taken into account. Firstly, 
we consider the hopping term, 



(ij)a 



(12) 



where t is the hopping parameter, (ij) means the nearest- 



neighbors, c' ia (ci a ) is the creation (annihilation) opera- 
tor of electrons at site i with spin a. With the basis 



3 



of {ai,a 2 }, we hereafter designate the coordinates of an 
unit cell as that of the black site at the left lower corner 
of the cell as shown in Fig. Q] The position of a site can 
then be denoted as (J, fj,) where j implies j-th unit cell 
and n = 1(2) corresponds to black (green) site. The elec- 
tron operator, for example the annihilation one, should 
be then denoted as Cj QjM . Expanding the operator with 
the plane waves, we have 



lv Cje 



exp(ifc • fj 



(13) 



with N as the number of total unit cells of the lattice, 
and fj as the position vector of the j-ih unit cell. Under 
such a convention, we can easily rewrite Hq in momentum 
space. The result is 




H 



k a 



(14) 



here the electron operators are given by spinors, tpt — 
(ct , ct ) with the first and second components denot- 

v ka,l' ka,2> 1 

ing electrons respectively at the black and green sublat- 
tices 
and 



tices, hz = £ k *i CJl t k2 tj2 w ith it's as the Pauli matrices, 



C k2 



-t(l + cos k x + cos k y ) 
-i(sin k x + sin k y ) 



(15) 
(16) 



with k x and k y as the components of the electron mo- 
mentum k in the basis of {bi,b2j-. 

The special feature of graphene is in the energy disper- 
sion. The eigenvalues of hj; are ie(fc) with 



e(k) 



£ kl 



el . 

k2 



(17) 



In Fig. [31 we show the energy dispersion of upper band. 
Clearly, e(fc) depends on k linearly only when k closes 
to ±(1,— l)27r/3 where it vanishes. In the continuous 
model, the energy dispersion is approximated by simple 
cones. 

We next consider the Coulomb interaction. The inter- 
action Vjj, u (Ri , Rj ) between two electrons respectively at 
positions (i, /i) and (j, v) is given by 



FIG. 3: (Color online) The energy of an electron in the upper 
band of non-interacting electrons as function of k 



L is the vector from the black site to green site in an 
unit cell as shown in Fig. 1. The on-site interaction 
U is the Coulomb repulsion between electrons of anti- 
parallel spins, leading to the short-range antiferromag- 
netic correlations (AFC). Since the AFC is not significant 
in graphene, U should not be too large. In our calcula- 
tion, we set U = 2e 2 /eL which is double of the nearest- 
neighbor interaction. For the long-range Coulomb inter- 
acting system, the final result should not sensitively de- 
pend on such a small but reasonable U. By taking into 
account only the charge fluctuations (with the spin fluc- 
tuations neglected), the interaction term of the Hamilto- 
nian in momentum space is given by 



Hi = 



-E 



(21) 



where nt 
. 9 



(nt^nt^) is the electron density operator, 
and iiq- is the Fourier component of the Coulomb interac- 
tion. Since the total charge of the electrons is neutralized 
by the background of positive charges, the q — term is 
excluded from the summation 
given in Appendix A. 



The elements of are 



y^{Ri,Rj) 
V 12 (Ri,Rj) 
V 21 (Ri,Rj) 



e\Ri-R 3 \ 

U, for R; = R 



e|i?i - Rj - L\ 



e\Ri — Ri 



for Ri Rj 



(18) 
(19) 
(20) 



where Raj) denotes the position of the i (j)-th unit cell, 
e is the static dielectric constant due to the screening 
by the electrons of carbon core and the substrate, and 



IV. GREEN'S FUNCTION 



The Green function of electrons is defined as 

G(k,T-r>) = -{Tt^Jt^IJt')), (22) 

where r is the imaginary time, and (T T ■ ■ ■) means the 
statistical average of the r-ordered product of the oper- 
ators. In the frequency space, G is given by 



G(k,iu>n) = [iu n - £,t ~ t,{k,iuj n )] \ 



(23) 
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FIG. 4: (Color online) The electron self-energy in the 
renormalized-ring-diagram approximation. The solid lines are 
Green's functions, the single and double dashed lines are re- 
spectively the bare and renormalized Coulomb interactions, 
and the bubble is the electron polarizability. 



where ^ = h^ — [i with \i the chemical potential, uj n is the 
fermionic Matsubara frequency, and S(fe, iw n ) is the self- 
energy. For brevity, we hereafter will use k = (fc, iui n ) for 
the arguments unless stated otherwise. Under the RRDA 
which is diagrammatically shown in Fig. 21 the elements 
of S(fc) are given by 

S^(fc) = -^E G ^( fc -9)^(9)> (24) 
i 

where [3 = 1/T with T as the temperature (in unit of 
fcs = 1) of the system, W^(q) is the element of the 
matrix of the screened interaction 

W(q) = [1 - t>( 9 )x(g)]- 1 F(g), (25) 

and the element of the polarizability x(q) is given by 

= ^ E G^(k + q)G^{k) (26) 

with q = (q, ifl m ) and Q m is the bosonic Matsubara fre- 
quency. The chemical potential /i is determined by the 
electron density (that is the electron number per site), 



n 



— ^Gn(fc)exp(iu n O+). (27) 



From equations l|23 p -(|27 p . the Green functions can 
be self-consistently determined. Using the Pade 
approximation,— we obtain the retarded self-energy 
Yt r (k,to) and then the Green function G r (fc, u>) under the 
analytical continuation iu> n — > u> + i0 + . 

In order to do an effective numerical calculation, it is 
necessary to get a clear understanding of the symmetry 
of the Green functions, and this is shown in Appendix B. 

On the other hand, we need to pay attention to the 
behavior of the screened interaction W(q). Since it ap- 
proaches to the bear Coulomb interaction v(q) in the 



limit iQ m — » oo, W{q) can be separated into two parts, 
W(q) = v(q) + Wn(q), where Wn{q) is the induced inter- 
action from the ring diagrams. The contribution of v(q) 
in the self-energy yields the exchange part. Here, an im- 
portant point is that the behavior of Wn(q) in the limit 
q — > by RRDA is very different from that by RPA. For 
the sake of illustration, we use the descriptions for v(q) 
and x(q) by their Pauli components. In the limit q — ► 0, 
we have 



v(q) 



r 



\0- + <ri) 



X.{l) -* Xo(0,iOm) +Xl(0,i^m) cr l 

where V — Aire 2 j y3ea. For Wn(q) in the same limit, we 
get 



W R (q) - - 



Q(Q 



-(l + tri) 



(28) 



with q m = -2T[ Xo (0,in m ) + Xi(0,iO m )]. In RPA, q m 
vanishes for m ^ 0. However, in RRDA, q m ^ 0, it means 
the dynamic screening effect in the long-wavelength in- 
teractions. 

As is well known, the ring diagrams in the RPA 
include the contribution of plasmon excitations. The 
plasmon frequency Qq oc ^/Q is determined by the 
long-wavelength behavior of the charge polarizabilty 

XRPA(q,i&m) = [Xo(<f,^m) + Xl (<f, i&m)]llPA OC Q 2 at 

q — > 0. This polarizability is calculated in the absence 
of interactions. Under the RRDA, the corresponding po- 
larization diagram of Xo(Qi^m) + Xi{Q> i^m) is calcu- 
lated with renormalized Green functions, and it does not 
have such a behavior as in RPA. The renormalized-ring- 
diagram summation does not result in the desired plas- 
mon excitations. The correct way to obtain the plasmon 
excitations is to calculate the two-particle Green func- 
tion in which the vortex corrections need to be consid- 
ered. Under the RRDA that is a conserving approxima- 
tion for the single-particle Green function, the kernel of 
the equation for the two-particle propagator is generated 
from the functional derivative of the self-energy diagrams 
with respect to the Green function.— The calculation of 
the two-particle propagator needs a more complicated 
mathematical procedure and is beyond the scope of the 
present approach. 



V. NUMERICAL METHOD 

Under the present approximation, the screening effect 
is negligible only at sufficient large Matsubara frequen- 
cies. This requires a considerable amount of numerical 
calculations. To save the computer time, the summa- 
tions over the Matsubara frequencies in Eqs. (|24[) and 
([26)1 need to be performed with special method. We 
have developed a super-high-efficiency algorithm for the 
series summations^ In the present calculation, we have 
used the parameters [h, L, M) — [2, 15, 5] for selection of 



5 



the Matsubara frequencies distributed in L successively 
connected blocks each of them containing M frequencies 
with h the integer-parameter that the stride in the ^-th 
block is h^~ x \ The total number of the frequencies se- 
lected here is L(M — 1) + 1 = 61. The largest number 
N c is about 2 L (M — 1) = 2 17 for the cutoff frequencies 
Q Nc = 2N c ttT and uj Nc = (2N C - i)wT. For the low- 
est temperature considered here, T/t = 0.01, we have 
k>N c /t ~ 8235. At the frequencies larger than the cutoff, 
Wn{q) is negligible small. For illustrating the numerical 
method, an example for calculation of the element Xwil) 
is given in Appendix C. 

On the other hand, the momentum integrals in Eqs. 
and (f2"6")) are convolutions. These integrals can be 
efficiently carried out by Fourier transforms. Here again, 
we should pay attention to Wn{q). It should be care- 
fully transformed from q-space to r-space since it has a 
sharp peak at q — as indicated by Eq. (f2"8")) . This 
long-wavelength singular part should be subtracted from 
Wn{q) and can be treated especially. It saves computer 
time to perform the transform of this singular part in 
the orthogonal coordinate system because where it is 
isotropic and the Q-integral within a circle close to the 
origin can be reduced to a one-dimensional one. Only 
within this circle, we need a very fine mesh for the in- 
tegral. In the rest part of the hexagon Brillouin zone, 
one can use a crude mesh doing the two-dimensional in- 
tegral because where the integrand is not singular. The 
remained part of Wu(q) should be a regular function ex- 
cept there may be some undulations close to q = due 
to the subtraction. 

Another point we should take care about is that the 
Green function varies drastically around and close to the 
Fermi surface at low temperatures. Since the band struc- 
ture is not flat at the Fermi energy, the mesh for q-space 
integral should be fine enough around and close to the 
Fermi surface. We here give an example for sampling 
the points in momentum space. This example is for the 
zero-doping. Under the same consideration, the sampling 
for finite doping cases can be planed similarly. We divide 
the range [0, n] in each axis in momentum space into four 
blocks shown as in Fig. [5l There are Ni equal meshes in 
the i-th block and iVj's are given by 

{N U N 2 ,N 3 ,N 4 } = {6,8,40,4}. 

The finest mesh is for the third block 22C] which is 
centered with 2ir/3. In the Brillouin zone, this leads to 
a very fine mesh around the Dirac points ±(1, — l)27r/3. 
The second finest mesh is for the first block, which is 
designed for dealing with the undulations of the remained 
part of Wn(q) (after the subtraction of the sharp peak) 
when it is Fourier transformed from g-space to r-space. 
The rest two blocks have relatively crude meshes because 
the Green function is smooth there. 

Our numerical algorithm considerably saves the com- 
puter memory and time. The similar method and its 
accuracy have been demonstrated by a recent study on 
the two-dimensional electron system with infinite band 



40 



FIG. 5: For zero-doping calculation, the range [0, n] is divided 
into 4 blocks with partition points 7r/12, 77r/12 and 3n/4. 
Each number in the third line is the number of meshes in 
the corresponding block. 



width and long-range Coulomb interaction^ With our 
numerical algorithm, we have solved the above equations 
by iteration. 



VI. RESULTS 

The system is characterized by the coupling constant 
g that is defined by the ratio between the overall-average 
interaction energy e 2 /ea and the hopping energy t, 



e 

eat 



(29) 



The parameters t — 2.82 eV and a — 2.4 A are known 
from the experimental observations.— By choosing e ~ 
4, we have g ~ 0.5. Therefore, the graphene is a 
moderately-coupled Coulomb system. 

Firstly, we present the result for the spectral density 
that is defined by 



A(k,E) 



-ImTrG r (k,E) 
'-lmG r n (k,E), 



(30) 



where G\ 1 {k,E) = G'2 2 (k,E). In the non-interacting 
case, A(k, E) reduces to the 5-functions representing the 
energy dispersions of the two bands. In the present case, 
the energy levels are broadened because of the many- 
body effect. Shown on the left panel in Fig. [B]is an inten- 
sity map of the spectral density in the energy-momentum 
plane at the doping concentration c = n — 1 = 0.02 and 
temperature T/t = 0.02. This map exhibits the energy 
distribution of the states as a function of momentum 
along the high symmetry directions T-M-K-T in the Bril- 
louin zone. For comparison, the free-particle energy dis- 
persion — e(k) — n° (with /z° the chemical potential of 
the non-interacting system) is also depicted as the solid 
curve. Clearly, because of the Coulomb effect, the en- 
ergy distribution has finite width, implying the finite life 
times of quasiparticlc. In addition, the scale of the en- 
ergy band is enlarged. The dashed curve is a rescale of 
the solid one. With comparing to the non-interacting dis- 
persion, the energy band is magnified with a factor about 
1.1. Actually, the exchange self-energy results in an ad- 
ditional hopping of the electrons, which renormalizes the 
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FIG. 6: (Color online) Left panel: Spectral density as a func- 
tion of energy along the high symmetry directions in the Bril- 
louin zone. Right panel: Spectral density at zero energy show- 
ing the Fermi surface structure. The doping concentration of 
the electrons is c = 0.02, and the temperature is T/t = 0.02. 

kinetic energy^ To see this, we express the exchange 
self-energy in real space, 

Ef 3 .=-(ct aCfa K (31) 

where Vij is the Coulomb interaction between electrons 
at sites i and j. [Here i and j mean the positions of 
the lattice sites.] With the quantity X?-, one can de- 
fine the additional hopping parameter, tfj = — £?■. The 
most sizable additional hopping is the one between the 
nearest-neighbor (NN). At c = 0.02, and g = 0.5, we 
find t x NN = 0.224, and the magnitude of the next NN 
one is two order less than t x NN . At large distance, tfj is 
very small. As a total effect, such the additional hop- 
ping terms enlarge the original energy band. On the 
right panel of Fig. [51 we show the map of the spectral 
density obtained by integration over the energy window 
of 0.064 around the chemical potential. The orbits of 
strong intensity correspond to the Fermi surfaces which 
are apparent not circles as compared with those from the 
simplified Dirac cone model. The structure of the Fermi 
surfaces is symmetrical under any rotation of angle 7r/3 
around the origin. All these results are comparable with 
the ARPES experimental observations i 15 ' 17 For c=0.0, 
we expect that the Fermi surfaces in Fig. [5] shrink into 
Dirac points. 

The broadening of energy distribution of quasipar- 
ticle is described by the imaginary part of the self- 
energy. To be specific, we analyze the Green function 
G^(k, E). This function can be divided into two parts, 
[Gl(k, E) + G\{k, E)]/2 with 



S(k, E) = {[ eg 12 + YT 12 (%, E)][e % 2l + YT 2l {k, E)} 1 ' 2 

where e^ 12 = ei n = e R>1 -ie^ 2 , £g(fc, E) = E^(fc, E) = 

Yi22{k,E) : and E^(fc, E)'s are the components of the 
retarded self-energy matrix. The two Green functions 



c = 0.02, g = 0.5, T/t = 0.02 
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FIG. 7: (Color online) Imaginary part of the self-energy 
EJ^ i(k,E) as function of E at doping concentration c = 0.02 
and T/t = 0.02 for the upper and lower bands. The result 
from the RPA for the upper band is also plotted here for 
comparison. 

G u and Gi can be considered as for the electrons at the 
upper band and lower band, respectively. Corresponding 
to G r u ; (fc, E), we define the retarded self-energies for the 
electrons respectively at upper and lower bands as 

Z r U!l (k,E) = X r (k,E)±S(k,E). (32) 

In Fig. the imaginary parts of the self-energies 
; (k, E) are presented as functions of E at a Fermi 

momentum k w (— 0.87T, 0.67r) for the system at doping 
concentration c = 0.02 and the temperature T/t = 0.02. 
At small energy, Im£^ (£;,£?) is a quadratic function of 
E. The value at E = is very small and should van- 
ish at zero temperature. At a region of larger energy, 
ImI7 u (k,E) seems to be linearly depended on E. The 
magnitude of Im£y(fc,_E) is, however, small compared 
with ImE^(fc, E). The reason is clear. At the Fermi sur- 
face, the states of energy E is far from the correspond- 
ing states of the same momentum k at the lower band. 
Therefore, the self-energy T,J(k,E) is small. In Fig. [71 
the RPA result for the self-energy of upper band is also 
shown for comparison. The magnitude of Im£^(fc, E) by 
RPA is larger than that by the RRDA. But, the quadratic 
-E-dependence of ImS^(fc, E) at small E is also reflected 
by the RPA calculation. 

In Fig. the results for Im£^(fc, E) are depicted 
at low temperatures and at zero doping concentration. 
Since each Fermi surface is a Dirac point in this case, 
the Fermi momentum for Fig. [5] is chosen as k = 
(— 27t/3, 27t/3). Firstly, at the Dirac point, because the 
state is the common state of the upper and lower bands, 
the self-energies of both bands coincide. At small en- 
ergy, apart from a small value at E = due to the finite 
temperature effect, all the results for ImS^(fc, E) show a 
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FIG. 8: (Color online) Imaginary part of the self-energy 
E) as function of E at zero doping concentration for 
several different temperatures. 



FIG. 9: (Color online) RPA result for the imaginary part 
of the self-energy Yi T u l (k,E) as function of E at zero doping 
concentration. 



quadratic dependence on E. Out of the small energy re- 
gion, ImE^(fc, E) linearly depends on E to certain limit. 
For the purpose of comparison, The RPA results obtained 
from our approach at the zero doping concentration are 
exhibited in Fig. M In the limit T -> 0, lmE r u l (k, E) 
appears to be linearly depended on \E\. This is in agree- 
ment with the analytical result obtained from the contin- 
uous model based on the RPA. 6 Therefore, the system at 
zero doping is referred^ to behave like a marginal Fermi 
liquid 18 with the imaginary part of the self-energy going 
linear in energy near the chemical potential. 

However, our numerical results based on the RRDA 
demonstrate that the imaginary parts of the self-energies 
fixed at the Fermi momentum are always varying as 
quadratic in energy close to the Fermi level, regardless 
the system is doped or not doped. This feature indicates 
that the quasiparticle in graphene behave like a moder- 
ately correlated Fermi liquid. 

With the retarded Green function, we can also calcu- 
late the density of states (DOS) defined as 

P( E )^-^J2 lmG ^ E ^ (33) 

k 

Shown in Fig. |TD] is the result for p(E) at zero doping 
concentration and T/t = 0.02. The non-interacting coun- 
terpart p°(E) is also depicted for comparison. Since the 
energy linearly depends on the magnitude of the momen- 
tum near the Dirac points, p (E) is proportional to E at 
small E. The two peaks come from the van Hove sin- 
gularity because the energy bands are flat at E — ±t as 
shown in Fig. [3l Under the Coulomb interaction, the 
spectral density of quasiparticle is broadened. This re- 
sults in lowering the density of states and smearing the 
peaks. With comparing to p°(E), the two peaks in p(E) 
shift to larger energy because of the energy band enlarged 
by the exchange interaction. 
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FIG. 10: (Color online) Density of states p(E) as function of 
E at zero doping concentration. p°(E) is the non- interacting 
counterpart. 

We next consider the distribution function n(k). For 
the graphene, n(k) should be defined as 

n[k) = \Ys Tr( ^' exp(iw„0 + ). (34) 

n 

As we have encountered in Eq. (|30jl . n(k) is determined 
only by Gn(fc, iu> n ). Further more, since the Green func- 
tion Gii(fc,za>„) can be divided into parts of the up- 
per and lower bands, we therefore define the distribution 
functions for these two bands as 

n Ul i(k) = ^^G n j(fc,ia;„)exp(zw n + ). (35) 

n 

The total distribution function is given by n(k) — n u (k) + 
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FIG. 11: (Color online) Distribution functions at doping con- 
centration c = 0.02 and T/t = 0.02. 



ni(k). In Fig. [TTJ we exhibit the result of the distribu- 
tion functions n u ,i{k) at temperature T/t = 0.02 and 
doping concentration c = 0.02. In the inner Fermi area, 
the upper band is almost fully occupied with the elec- 
trons. The distribution n u (k) drops drastically at the 
Fermi surface. The occupation at the lower band is not 
flat. Close to each Dirac point (the corner of the hexagon 
Brillouin zone), the behavior of ni(k) looks like a sink. 
The depression of n/(fc) comes from two aspects. One 
is the temperature effect. The distribution ni(k) reaches 
its minimum at the Dirac points because where the en- 
ergy (Dirac energy) is the highest for states in the lower 
band. Since the doping concentration is small, the Dirac 
energy is close to the Fermi level. Therefore, ni(k) has 
an apparent drop at the Dirac points. Another reason is 
due to the many-body effect. Because of the Coulomb 
interaction, an amount of electrons can be redistributed 
from the lower band to upper band. 

At zero doping, each Fermi surface shrinks to a point. 
The distribution functions are shown in Fig. [12] The dis- 
tribution n u {k) at upper band concentrates at the Dirac 
points. At lower band, there is a obvious depression of 
ni{k) at the point. To see the many-body and temper- 
ature effects in the zero-doping case, we show the to- 
tal distribution n(k) at and close to the Dirac points as 
functions of temperature. Very close to the Dirac points, 
n{k) increases drastically with decreasing temperature. 
At zero temperature, the total distribution at the Dirac 
points 1 < n(fco) < 2 can be expected. This is very differ- 




FIG. 12: (Color online) Distribution functions at zero doping 
concentration and T/t = 0.02. 



ent from the non-interacting distribution. The latter is 
constant 1 because the zero-energy levels of both bands 
are half occupied and at the momentum different from 
the Dirac points only the lower band is fully occupied. 
Under the Coulomb interactions, some of the electrons 
around each Dirac point in the lower band are gathered 
to the upper band close to the Dirac point, resulting in 
a higher distribution at and close to the point. From the 
increasing tendencies of n(fco) and n(ko + A/co), we can 
infer there is an abrupt drop in n(k) close to each Dirac 
point. Therefore, the zero-doping distribution function 
at zero temperature is consistent with the Fermi liquid 
behavior. 

For the negative doping, the Fermi surface opens again. 
In this case, the Fermi level is at the lower band. Since 
the upper and lower bands are symmetric about the zero 
energy, at zero temperature, the electron distribution at 
the upper band corresponds to the hole distribution at 
the lower band, and vice versa. Therefore, the distribu- 
tions at negative doping can be obtained from that of 
positive doping. For example, at c = -0.02, by flipping 
Fig. [IT] upside down, we obtain the image of distribu- 
tions. The shape of each Fermi surface is the same as 
that at c = 0.02. However, the Fermi area is outside of 
the surface where the lower band is near fully occupied. 

Finally, we give the ground-state energy per electron 
e - At low temperatures, T/t < 0.1, the numerical results 
for the energy per electron are almost a constant. By 
extrapolation, we then obtain eo- The results for the two 
doping cases are, eo = — l.lOt for c = 0, and eq = — 1.092 
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FIG. 13: (Color online) Temperature dependence of the total 
distribution at and close to the Dirac points at zero doping 
concentration. 



for c = 0.02. 



VII. SUMMARY 

In summary, using the tight-binding model with long- 
range Coulomb interactions defined in a honeycomb lat- 
tice, we have presented the Green-function formulation 
for the electron in graphcnc. The interactions between 
electrons are treated with the renormalized-ring-diagram 
approximation. The integral equations for determining 
the Green function are solved self-consistently using our 
high-efficiency numerical algorithm. The obtained spec- 
tral densities are comparable with the experimental ob- 
servations. Since the imaginary part of the self-energy 
of the electron Green function fixed at the Fermi mo- 
mentum varies as quadratic in energy near the chemical 
potential for both doped and undoped systems, we con- 
clude that electrons in graphene follow the Fermi liquid 
like behavior. In addition, we also calculated the density 
of states and the distribution functions. 



Acknowledgments 

This work was supported by a grant from the Robert 
A. Welch Foundation under No. E-1146, the TCSUH, 
and the National Basic Research 973 Program of China 
under grant No. 2005CB623602. 

APPENDIX A: Fourier transform of Coulomb 
interaction between electrons on the honeycomb 
lattice 

In this appendix, we present the Fourier transform of 
the Coulomb interaction between electrons on the honey- 



comb lattice. Because it is long-range interaction, the ac- 
curate result cannot be obtained from the direct summa- 
tion by definition with numerical calculation. We must 
seek a fast converging scheme for summation over the 
lattice sites. To do this, we separate the interaction 



V(R) 



|r, for i?^0 
U, for R = 



(36) 



into short-range and long-range parts, V{R) — Vs(R) + 
V L (R) with 



V L (R) 



eVa 2 + R 2 ' 



(37) 



where a is a free parameter of positive quantity, and 
V S (R) = V(R) - V L (R). At R > a, V S (R) = 0(iT 3 ), 
therefore it is a short-range interaction. Its Fourier trans- 
form can be obtained by the direct summation, 



Vs(Q) = ^V s (R j )e- i & ii i 



(38) 



On the other hand, Vl(R) is long ranged. Direct summa- 
tion for the Fourier transform of it converges very slow. 
However, this function can be expressed as 



Vl(R) 



dQ 
(2*Y 



HQ) 



,iQR 



(39) 



where the Q-integral is over the whole momentum space, 
and 



HQ) 



27TC 2 



-aQ 



(40) 



The integral in Eq. (|39[) can be written in a summation 
over the integrals each of them over a Brillouin zone. 
Shifting all these Brillouin zones to the first Brillouin 
zone by the corresponding reciprocal lattice vectors Q„'s, 
we have 



v L (R) = J2 



BZ 



dQ 

(27T) 1 



H\Q + Qn\)< 



i(Q+Q n )-R 



(41) 



The order of integral and summation can be changed. 
Express R as R = Rj + Z where Rj is the position of 
j-th unit cell of the honeycomb lattice, and Z is a vector 
within the unit cell. We get 



VU\Rj 



BZ 



dQ 
(2n) 2 



Y,H\Q+Qn\)e l{Q+Q " yZ e lQ - R 

n 

(42) 

where use of exp(iQ„ • Rj) = 1 has been made. From this 
equation, we recognize the Fourier transform of the func- 
tion F{Rj 1 Z) = VtiRj + Z) defined on the honeycomb 
lattice with Z a parameter, 



F(Q, Z) = 1=Y,H\Q + Qn\)e 



i{Q+Q n )-Z 



(43) 
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where the factor 2>/3 comes from the ratio between the 
area of the first Brillouin zone of the honeycomb lattice 
and that of the square lattice (2tt) 2 . Since <fi(Q) decreases 
exponentially at large \Q n \, the n-summation converges 
very fast. 

We are now ready to express the elements of the in- 
teraction matrix v(q) = v s (q) + v L (q). With the above 
functions, wc have 



(44) 



v?2(?) = v^(q) = J2Vs(R j -L)e- i ^, (45) 



ViM = v£ 2 (q) = F(Mq,0), 
vUo) = v£(q)=F(Mq,-L), 



(46) 
(47) 



where Mq maps the vector q under the basis {bi,b2} 
into the one in the orthogonal coordinate system. 

APPENDIX B: Symmetry of the Green func- 
tion defined on the honeycomb lattice 

For doing numerical calculation, it is necessary to un- 
derstand the symmetry of the Green function. We start 
the discussion with the definition of the Green function 
in real space, 

G^{r i -r j ,r-T') = -{T T c ia ^T)c] atV {r')), (48) 

which describes a particle of spin-a propagating from po- 
sition (j, v) to (i, fi). Firstly, this function has the follow- 
ing property, 

G„„(n -fj,T- t') = G^ifj -n.T-r') (49) 

because that the configuration for a particle propagating 
from (j, v) to (i, fi) is the same as in the inverse process. 
Therefore, the diagonal Green functions are even under 
r — ► — f. Denoting Go = Gn = G22, we obtain 

G (r,r) =G (-r,T). (50) 

On the other hand, for the off-diagonal part, we get 

G 12 (f,T) = G 21 (-f,r). (51) 

Though the off-diagonal Green functions have no definite 
parity, they can be separated into even and odd functions. 
For G12 (r , r) , we have 



G 12 (f,r)=G 1 (f,r)-iG 2 (f,T) 



(52) 



where Gi (r, r) and G2 (r , t) are even and odd functions of 
f, respectively. Now, we can express the Green function 
matrix in terms of Pauli matrix, 

G(f,T)=G Q (f,T)ao+G 1 (f,T)a 1 +G 2 (r f ,T)a2 (53) 



where the Pauli components as functions of f have their 
definite parities. Obviously, in momentum space, as func- 
tions of k, they have the same parities. Using the prop- 
erty of parity, the Fourier expansions of these components 
are given by 



Go,i(r,r) 



G 2 (f,r) 



H 



dk 

—Go A (k,T)co S (k-r) (54) 



JH 



dk 
2^2 



G 2 (£:,T)sin(fc-r) (55) 



where the fc-integrals are over the half Brillouin zone: 
< k x < 7r and — ir < k y < ir. Further more, in the 
half Brillouin zone, we can separate the Green functions 
G a (k, t)'s into even and odd functions of k y , G*(fc, t) 
with the superscripts ± denoting the parities, 

G±(fc x ,fc y ,r) = [G a {k x ,k y ,r)±G a {k x ,-k v ,T)]/2. (56) 

Correspondingly, in real space, G CT (r, r)'s are separated 
into even and odd functions of r y . We have 

f dk 

G oi(r,T) = / — G f 1 (fc,T)cos(fc x r a; )cos(fcyrj / ) 

Jl 71" 

f dk 

G o,i(^ T ) = / -^G^ 1 (k,T)sm(k x r x )sm(k v r y ) 
f dk 

G^{r,r) = i I —~G%(k,T)sm(k x r x )cos(k y r y ) 

Jl 7T Z 



G 2 {f 



f dk 



G 2 (fc, t) cos(k x r x ) sm(k y r y ) 



where the /c-integrals are now over the first quadrant of 
the first Brillouin zone: < k x < tt and < k y < it. 

We next consider the property of the Green functions 
under the exchange of the coordinates (r Xl r y ) — > (r yi r x ). 
Since the system is symmetric under this exchange, we 
have 



G in/ {t x , T y , t) — G (r y , v x , t) . 



(57) 



From the above definitions, one can obtain that the Pauli 
components G^ are symmetric about the exchange but 

G 2 (r x ,r y ,r) = G2(r y ,r x ,T) 
G2(r x ,r y ,r) = G^(r y ,r x ,T). 

Making use of these symmetries in a numerical process, 
we need to calculate these Green functions only in half 
of the first quadrant in both real and momentum spaces. 

Finally, the symmetry related to the time reversion is 

Gt*(k,iu n ) =Gt(k,-iu n ) for a = 0,1, 2 (58) 

which is obtained from the definition. 

The above discussion of symmetries for the Green func- 
tions applies to the Coulomb interactions, polarizabili- 
ties, and the self-energies as well. 
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APPENDIX C: An example for calculation of 

For the readers' convenience, we here give an example 
for calculation of X^il)- F° r the briefness, we present 
only the result for X^io) = Xo(q)- The results for other 
elements and the method for calculation of the self-energy 
elements E A11/ (fc) can be obtained immediately with the 
same consideration. As has been mentioned in the main 
text, the momentum convolution of two Green functions 
in Xo(?) can be evaluated by Fourier transform. In real 
space, it is given by 

2 

Xo(f, iCl m ) = — 22 Go(r, iuJ n + i£l m )G (r, iu> n ) (59) 

" n 

where the n-summation for the Matsubara frequencies is 
over (— oo, oo). Using the property of the Green function, 
Gg(r, iuj n ) = Go(r, —iuj n ), we take the summation only 
for n > 0, 

. oo 

4 x 

Xo(r,iO, m ) = -Re{ ) j G {r, iuj n + iQ m )G (r, iu n ) 

P 71=1 

[m/2] 

+ 2J G Q {f,in m - ico n )G (f, -iu n ) 

71=1 

+ \G (r,iu n )\y2\ri^Z 2 }. 



where [m/2] is the integer part of m/2. Applying our 
rule for the series summation, we get 

Xo{f, itt m ) = -Re{ 2 J WpG (r, iuj p + iQ m )G (f, imp) 
p 

+ w l p m/2] G (r, ifl m - iuj p )G Q (r, -iuj p ) 
p 

+ \G a (r^)\y2\rL m :Z 2 } 
+ $Xo(r,itt m ) 

where the summation is over the selected points {p} with 

w p and w p " 1 ^ 2 ^ the corresponding weights^ and the last 
term is the contribution from terms beyond the cutoff N c 
for the Matsubara frequency, 

{N c +m 
-^r E a^i ifm>0 
n=l 

(60) 

The present expression for vn is essentially the same as 
that in Appendix B of Ref . [lj except an additional factor 
2 stemming from the degree of spin freedom in the present 
case and a misprinting in Ref. Il4l . The notation of taking 
real part of the result should be assigned in the expression 
of x m Appendix B of Ref. [13 . . 
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